/******************************************************************************
 *
 * MantaFlow fluid solver framework
 * Copyright 2011 Tobias Pfaff, Nils Thuerey 
 *
 * This program is free software, distributed under the terms of the
 * Apache License, Version 2.0 
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Meshes
 *
 * note: this is only a temporary solution, details are bound to change
 *        long term goal is integration with Split&Merge code by Wojtan et al.
 * 
 ******************************************************************************/

#include "mesh.h"
#include "integrator.h"
#include "mantaio.h"
#include "kernel.h"
#include "shapes.h"
#include "noisefield.h"
//#include "grid.h"
#include <stack>
#include <cstring>

using namespace std;
namespace Manta {

Mesh::Mesh(FluidSolver* parent) : PbClass(parent) {  
}

Mesh::~Mesh()
{
	for(IndexInt i=0; i<(IndexInt)mMeshData.size(); ++i)
		mMeshData[i]->setMesh(NULL);

	if(mFreeMdata) {
		for(IndexInt i=0; i<(IndexInt)mMeshData.size(); ++i)
			delete mMeshData[i];
	}
}

Mesh* Mesh::clone() {
	Mesh* nm = new Mesh(mParent);
	*nm = *this;
	nm->setName(getName());
	return nm;
}

void Mesh::deregister(MeshDataBase* mdata) {
	bool done = false;
	// remove pointer from mesh data list
	for(IndexInt i=0; i<(IndexInt)mMeshData.size(); ++i) {
		if(mMeshData[i] == mdata) {
			if(i<(IndexInt)mMeshData.size()-1)
				mMeshData[i] = mMeshData[mMeshData.size()-1];
			mMeshData.pop_back();
			done = true;
		}
	}
	if(!done)
		errMsg("Invalid pointer given, not registered!");
}

// create and attach a new mdata field to this mesh
PbClass* Mesh::create(PbType t, PbTypeVec T, const string& name) {
#	if NOPYTHON!=1
	_args.add("nocheck",true);
	if (t.str() == "")
		errMsg("Specify mesh data type to create");
	//debMsg( "Mdata creating '"<< t.str <<" with size "<< this->getSizeSlow(), 5 );

	PbClass* pyObj = PbClass::createPyObject(t.str() + T.str(), name, _args, this->getParent() );

	MeshDataBase* mdata = dynamic_cast<MeshDataBase*>(pyObj);
	if(!mdata) {
		errMsg("Unable to get mesh data pointer from newly created object. Only create MeshData type with a Mesh.creat() call, eg, MdataReal, MdataVec3 etc.");
		delete pyObj;
		return NULL;
	} else {
		this->registerMdata(mdata);
	}

	// directly init size of new mdata field:
	mdata->resize( this->getSizeSlow() );
#	else
	PbClass* pyObj = NULL;
#	endif
	return pyObj;
}

void Mesh::registerMdata(MeshDataBase* mdata) {
	mdata->setMesh(this);
	mMeshData.push_back(mdata);

	if( mdata->getType() == MeshDataBase::TypeReal ) {
		MeshDataImpl<Real>* pd = dynamic_cast< MeshDataImpl<Real>* >(mdata);
		if(!pd) errMsg("Invalid mdata object posing as real!");
		this->registerMdataReal(pd);
	}
	else if( mdata->getType() == MeshDataBase::TypeInt ) {
		MeshDataImpl<int>* pd = dynamic_cast< MeshDataImpl<int>* >(mdata);
		if(!pd) errMsg("Invalid mdata object posing as int!");
		this->registerMdataInt(pd);
	}
	else if( mdata->getType() == MeshDataBase::TypeVec3 ) {
		MeshDataImpl<Vec3>* pd = dynamic_cast< MeshDataImpl<Vec3>* >(mdata);
		if(!pd) errMsg("Invalid mdata object posing as vec3!");
		this->registerMdataVec3(pd);
	}
}
void Mesh::registerMdataReal(MeshDataImpl<Real>* pd) { mMdataReal.push_back(pd); }
void Mesh::registerMdataVec3(MeshDataImpl<Vec3>* pd) { mMdataVec3.push_back(pd); }
void Mesh::registerMdataInt (MeshDataImpl<int >* pd) { mMdataInt .push_back(pd); }

void Mesh::addAllMdata() {
	for(IndexInt i=0; i<(IndexInt)mMeshData.size(); ++i) {
		mMeshData[i]->addEntry();
	}
}

Real Mesh::computeCenterOfMass(Vec3& cm) const {
	
	// use double precision for summation, otherwise too much error accumulation
	double vol=0;
	Vector3D<double> cmd(0.0);
	for(size_t tri=0; tri < mTris.size(); tri++) {
		Vector3D<double> p1(toVec3d(getNode(tri,0)));
		Vector3D<double> p2(toVec3d(getNode(tri,1)));
		Vector3D<double> p3(toVec3d(getNode(tri,2)));
		
		double cvol = dot(cross(p1,p2),p3) / 6.0;        
		cmd += (p1+p2+p3) * (cvol/4.0);
		vol += cvol;
	}
	if (vol != 0.0) cmd /= vol;    
	
	cm = toVec3(cmd);     
	return (Real) vol;
}

void Mesh::clear() {
	mNodes.clear();
	mTris.clear();
	mCorners.clear();    
	m1RingLookup.clear();
	for(size_t i=0; i<mNodeChannels.size(); i++)
		mNodeChannels[i]->resize(0);
	for(size_t i=0; i<mTriChannels.size(); i++)
		mTriChannels[i]->resize(0);

	// clear mdata fields as well
	for (size_t i=0; i<mMdataReal.size(); i++)
		mMdataReal[i]->resize(0);
	for (size_t i=0; i<mMdataVec3.size(); i++)
		mMdataVec3[i]->resize(0);
	for (size_t i=0; i<mMdataInt.size(); i++)
		mMdataInt[i]->resize(0);
}

Mesh& Mesh::operator=(const Mesh& o) {
	// wipe current data
	clear();
	if (mNodeChannels.size() != o.mNodeChannels.size() ||
		mTriChannels.size() != o.mTriChannels.size())
		errMsg("can't copy mesh, channels not identical");
	mNodeChannels.clear();
	mTriChannels.clear();
	
	// copy corner, nodes, tris
	mCorners = o.mCorners;
	mNodes = o.mNodes;
	mTris = o.mTris;
	m1RingLookup = o.m1RingLookup;
	
	// copy channels
	for(size_t i=0; i<mNodeChannels.size(); i++)
		mNodeChannels[i] = o.mNodeChannels[i];
	for(size_t i=0; i<o.mTriChannels.size(); i++)
		mTriChannels[i] = o.mTriChannels[i];
	
	return *this;
}

void Mesh::load(string name, bool append) {
	if (name.find_last_of('.') == string::npos)
		errMsg("file '" + name + "' does not have an extension");
	string ext = name.substr(name.find_last_of('.'));
	if (ext == ".gz") // assume bobj gz
		readBobjFile(name, this, append);
	else if (ext == ".obj")
		readObjFile(name, this, append);
	else
		errMsg("file '" + name +"' filetype not supported");


	// dont always rebuild...
	//rebuildCorners();
	//rebuildLookup();
}

void Mesh::save(string name) {
	if (name.find_last_of('.') == string::npos)
		errMsg("file '" + name + "' does not have an extension");
	string ext = name.substr(name.find_last_of('.'));
	if (ext == ".obj")
		writeObjFile(name, this);
	else if (ext == ".gz")
		writeBobjFile(name, this);
	else
		errMsg("file '" + name +"' filetype not supported");
}

void Mesh::fromShape(Shape& shape, bool append) {
	if (!append)
		clear();
	shape.generateMesh(this);
}

void Mesh::resizeTris(int numTris) {
	mTris .resize(numTris );
	rebuildChannels();
}
void Mesh::resizeNodes(int numNodes) {
	mNodes.resize(numNodes);
	rebuildChannels();
}

//! do a quick check whether a rebuild is necessary, and if yes do rebuild
void Mesh::rebuildQuickCheck() {
	if(mCorners.size() != 3*mTris.size())
		rebuildCorners();
	if(m1RingLookup.size() != mNodes.size())
		rebuildLookup();
}

void Mesh::rebuildCorners(int from, int to) {
	mCorners.resize(3*mTris.size());
	if (to < 0) to = mTris.size();        
	
	// fill in basic info
	for (int tri=from; tri<to; tri++) {
		for (int c=0; c<3; c++) {            
			const int idx = tri*3+c;
			mCorners[idx].tri = tri;
			mCorners[idx].node = mTris[tri].c[c];
			mCorners[idx].next = 3*tri+((c+1)%3);
			mCorners[idx].prev = 3*tri+((c+2)%3);
			mCorners[idx].opposite = -1;
		}
	}
	
	// set opposite info
	int maxc = to*3;
	for (int c=from*3; c<maxc; c++) {  
		int next = mCorners[mCorners[c].next].node;
		int prev = mCorners[mCorners[c].prev].node;
		
		// find corner with same next/prev nodes
		for (int c2=c+1; c2<maxc; c2++) {
			int next2 = mCorners[mCorners[c2].next].node;
			if (next2 != next && next2 != prev) continue;
			int prev2 = mCorners[mCorners[c2].prev].node;
			if (prev2 != next && prev2 != prev) continue;
			
			// found
			mCorners[c].opposite = c2;
			mCorners[c2].opposite = c;
			break;
		}
		if (mCorners[c].opposite < 0) {
			// didn't find opposite
			errMsg("can't rebuild corners, index without an opposite");
		}
	}    
	
	rebuildChannels();
}

void Mesh::rebuildLookup(int from, int to) {
	if (from==0 && to<0) m1RingLookup.clear();
	m1RingLookup.resize(mNodes.size());
	if (to<0) to = mTris.size();
	from *=3; to *= 3;
	for (int i=from; i< to; i++) {
		const int node = mCorners[i].node;
		m1RingLookup[node].nodes.insert(mCorners[mCorners[i].next].node);
		m1RingLookup[node].nodes.insert(mCorners[mCorners[i].prev].node);
		m1RingLookup[node].tris.insert(mCorners[i].tri);
	}
}

void Mesh::rebuildChannels() {
	for(size_t i=0; i<mTriChannels.size(); i++)
		mTriChannels[i]->resize(mTris.size());
	for(size_t i=0; i<mNodeChannels.size(); i++)
		mNodeChannels[i]->resize(mNodes.size());   
}

KERNEL(pts) returns(vector<Vec3> u(size))
vector<Vec3> KnAdvectMeshInGrid(vector<Node>& nodes, const FlagGrid& flags, const MACGrid& vel, const Real dt) {
	if (nodes[idx].flags & Mesh::NfFixed) 
		u[idx] = 0.0;
	else if (!flags.isInBounds(nodes[idx].pos,1)) 
		u[idx] = 0.0;
	else 
		u[idx] = vel.getInterpolated(nodes[idx].pos) * dt;
}

// advection plugin
void Mesh::advectInGrid(FlagGrid& flags, MACGrid& vel, int integrationMode) {
	KnAdvectMeshInGrid kernel(mNodes, flags, vel, getParent()->getDt());
	integratePointSet( kernel, integrationMode);    
}

void Mesh::scale(Vec3 s) {
	for (size_t i=0; i<mNodes.size(); i++)
		mNodes[i].pos *= s;
}

void Mesh::offset(Vec3 o) {
	for (size_t i=0; i<mNodes.size(); i++)
		mNodes[i].pos += o;
}

void Mesh::rotate(Vec3 thetas) {
	// rotation thetas are in radians (e.g. pi is equal to 180 degrees)
	auto rotate = [&](Real theta, unsigned int first_axis, unsigned int second_axis) 
	{
		if (theta == 0.0f)
			return;

		Real sin_t = sin(theta);
		Real cos_t = cos(theta);

		Real sin_sign = first_axis == 0u && second_axis == 2u ? -1.0f : 1.0f;
		sin_t *= sin_sign;

		size_t length = mNodes.size();
		for (size_t n = 0; n < length; ++n)
		{
			Vec3& node = mNodes[n].pos;
			Real first_axis_val = node[first_axis];
			Real second_axis_val = node[second_axis];
			node[first_axis]  = first_axis_val  * cos_t - second_axis_val * sin_t;
			node[second_axis] = second_axis_val * cos_t + first_axis_val  * sin_t;
		}
	};

	// rotate x
	rotate(thetas[0], 1u, 2u);
	// rotate y
	rotate(thetas[1], 0u, 2u);
	// rotate z
	rotate(thetas[2], 0u, 1u);
}

void Mesh::computeVelocity(Mesh& oldMesh, MACGrid& vel) {
	// Early return if sizes do not match
	if(oldMesh.mNodes.size() != mNodes.size())
		return;

	// temp grid
	Grid<Vec3> veloMeanCounter(getParent());
	veloMeanCounter.setConst(0.0f);

	bool bIs2D = getParent()->is2D();

	// calculate velocities from previous to current frame (per vertex)
	for (size_t i=0; i<mNodes.size(); ++i)
	{
		// skip vertices that are not needed for 2D 
		if(bIs2D && (mNodes[i].pos.z < -0.5f || mNodes[i].pos.z > 0.5f))
			continue;

		Vec3 velo = mNodes[i].pos - oldMesh.mNodes[i].pos;
		vel.setInterpolated(mNodes[i].pos, velo, &(veloMeanCounter[0]));
	}

	// discretize the vertex velocities by averaging them on the grid
	vel.safeDivide(veloMeanCounter);
}

void Mesh::removeTri(int tri) {
	// delete triangles by overwriting them with elements from the end of the array.
	if(tri!=(int)mTris.size()-1) {
		// if this is the last element, and it is marked for deletion,
		// don't waste cycles transfering data to itself,
		// and DEFINITELY don't transfer .opposite data to other, untainted triangles.

		// old corners hold indices on the end of the corners array
		// new corners holds indices in the new spot in the middle of the array
		Corner* oldcorners[3];
		Corner* newcorners[3];
		int oldtri = mTris.size()-1;
		for (int c=0; c<3; c++) {
			oldcorners[c] = &corners(oldtri,c);
			newcorners[c] = &corners(tri, c);
		}
		
		// move the position of the triangle
		mTris[tri] = mTris[oldtri];
		
		// 1) update c.node, c.opposite (c.next and c.prev should be fine as they are)
		for (int c=0; c<3; c++) {
			newcorners[c]->node = mTris[tri].c[c];
			newcorners[c]->opposite = oldcorners[c]->opposite;
		}
		
		//  2) c.opposite.opposite = c
		for (int c=0; c<3; c++) {
			if (newcorners[c]->opposite>=0)
				mCorners[newcorners[c]->opposite].opposite = 3*tri+c;        
		}
		
		// update tri lookup
		for (int c=0; c<3; c++) {
			int node = mTris[tri].c[c];
			m1RingLookup[node].tris.erase(oldtri);
			m1RingLookup[node].tris.insert(tri);
		}
	} 

	// transfer tri props
	for(size_t p=0; p < mTriChannels.size(); p++)
		mTriChannels[p]->remove(tri);
	
	// pop the triangle and corners out of the vector
	mTris.pop_back();
	mCorners.resize(mTris.size()*3);
}

void Mesh::removeNodes(const vector<int>& deletedNodes) {
	// After we delete the nodes that are marked for removal,
	// the size of mNodes will be the current size - the size of the deleted array.
	// We are going to move the elements at the end of the array
	// (everything with an index >= newsize)
	// to the deleted spots.
	// We have to map all references to the last few nodes to their new locations.
	int newsize = (int)(mNodes.size() - deletedNodes.size());
	
	vector<int> new_index (deletedNodes.size());
	int di,ni;
	for(ni=0; ni<(int)new_index.size(); ni++)
		new_index[ni] = 0;
	for(di=0; di<(int)deletedNodes.size(); di++) {
		if(deletedNodes[di] >= newsize)
			new_index[deletedNodes[di]-newsize] = -1;   // tag this node as invalid
	}
	for(di=0,ni=0; ni<(int)new_index.size(); ni++,di++) {
		// we need to find a valid node to move
		// we marked invalid nodes in the earlier loop with a (-1),
		// so pick anything but those
		while(ni<(int)new_index.size() && new_index[ni]==-1)
			ni++;

		if(ni>=(int)new_index.size())
			break;

		// next we need to find a valid spot to move the node to.
		// we iterate through deleted[] until we find a valid spot
		while(di<(int)new_index.size() && deletedNodes[di]>=newsize)
			di++;

		// now we assign the valid node to the valid spot
		new_index[ni] = deletedNodes[di];
	}
	
	// Now we have a map of valid indices.
	// we move node[newsize+i] to location new_index[i].
	// We ignore the nodes with a -1 index, because they should not be moved.
	for(int i=0; i<(int)new_index.size(); i++) {
		if(new_index[i]!=-1)
			mNodes[ new_index[i] ] = mNodes[ newsize+i ];
	}
	mNodes.resize(newsize);
	
	// handle vertex properties
	for (size_t i=0; i<mNodeChannels.size(); i++)
		mNodeChannels[i]->renumber(new_index, newsize);
		
	// finally, we reconnect everything that used to point to this vertex.
	for(size_t tri=0, n=0; tri<mTris.size(); tri++) {
		for (int c=0; c<3; c++,n++) {
			if (mCorners[n].node >= newsize) {
				int newindex = new_index[mCorners[n].node - newsize];
				mCorners[n].node = newindex;
				mTris[mCorners[n].tri].c[c] = newindex;
			}
		}
	}    
	
	// renumber 1-ring
	for(int i=0; i<(int)new_index.size(); i++) {
		if(new_index[i]!=-1) {
			m1RingLookup[new_index[i]].nodes.swap(m1RingLookup[newsize+i].nodes);
			m1RingLookup[new_index[i]].tris.swap(m1RingLookup[newsize+i].tris);
		}
	}    
	m1RingLookup.resize(newsize);
	vector<int> reStack(new_index.size());
	for(int i=0; i<newsize; i++) {
		set<int>& cs = m1RingLookup[i].nodes;
		int reNum = 0;
		// find all nodes > newsize
		set<int>::reverse_iterator itend = cs.rend();
		for (set<int>::reverse_iterator it = cs.rbegin(); it != itend; ++it) {
			if (*it < newsize) break;
			reStack[reNum++] = *it;
		}
		// kill them and insert shifted values
		if (reNum > 0) {
			cs.erase(cs.find(reStack[reNum-1]), cs.end());        
			for (int j=0; j<reNum; j++) {
				cs.insert(new_index[reStack[j]-newsize]);
#ifdef DEBUG
				 if (new_index[reStack[j]-newsize] == -1)
					errMsg("invalid node present in 1-ring set");
#endif
			}
		}
	}
}

void Mesh::mergeNode(int node, int delnode) {
	set<int>& ring = m1RingLookup[delnode].nodes;
	for(set<int>::iterator it = ring.begin(); it != ring.end(); ++it) {
		m1RingLookup[*it].nodes.erase(delnode);
		if (*it != node) {
			m1RingLookup[*it].nodes.insert(node);
			m1RingLookup[node].nodes.insert(*it);
		}
	}
	set<int>& ringt = m1RingLookup[delnode].tris;
	for(set<int>::iterator it = ringt.begin(); it != ringt.end(); ++it) {
		const int t = *it;
		for (int c=0; c<3; c++) {
			if (mCorners[3*t+c].node == delnode) {
				mCorners[3*t+c].node = node;
				mTris[t].c[c] = node;
			}
		}
		m1RingLookup[node].tris.insert(t);
	}
	for(size_t i=0; i<mNodeChannels.size(); i++) { 
		// weight is fixed to 1/2 for now
		mNodeChannels[i]->mergeWith(node, delnode, 0.5);
	}
}

void Mesh::removeTriFromLookup(int tri) {
	for(int c=0; c<3; c++) {
		int node = mTris[tri].c[c];
		m1RingLookup[node].tris.erase(tri);
	}
}

void Mesh::addCorner(Corner a) {
	mCorners.push_back(a);    
}

int Mesh::addTri(Triangle a) {
	mTris.push_back(a);
	for (int c=0;c<3;c++) {
		int node = a.c[c];
		int nextnode = a.c[(c+1)%3];
		if ((int)m1RingLookup.size() <= node) m1RingLookup.resize(node+1);
		if ((int)m1RingLookup.size() <= nextnode) m1RingLookup.resize(nextnode+1);
		m1RingLookup[node].nodes.insert(nextnode);
		m1RingLookup[nextnode].nodes.insert(node);
		m1RingLookup[node].tris.insert(mTris.size()-1);
	}
	return mTris.size()-1;
}

int Mesh::addNode(Node a) {
	mNodes.push_back(a);
	if (m1RingLookup.size() < mNodes.size())
		m1RingLookup.resize(mNodes.size());

	// if mdata exists, add zero init for every node
	addAllMdata();

	return mNodes.size()-1;
}

void Mesh::computeVertexNormals() {
	for (size_t i=0; i<mNodes.size(); i++) {
		mNodes[i].normal = 0.0;
	}
	for (size_t t=0; t<mTris.size(); t++) {
		Vec3 p0 = getNode(t,0), p1 = getNode(t,1), p2 = getNode(t,2);        
		Vec3 n0 = p0-p1, n1 = p1-p2, n2 = p2-p0;
		Real l0 = normSquare(n0), l1 = normSquare(n1), l2 = normSquare(n2);
		
		Vec3 nm = cross(n0,n1);
		
		mNodes[mTris[t].c[0]].normal += nm * (1.0 / (l0*l2));
		mNodes[mTris[t].c[1]].normal += nm * (1.0 / (l0*l1));
		mNodes[mTris[t].c[2]].normal += nm * (1.0 / (l1*l2));
	}
	for (size_t i=0; i<mNodes.size(); i++) {
		normalize(mNodes[i].normal);
	}
}

void Mesh::fastNodeLookupRebuild(int corner) {    
	int node = mCorners[corner].node;
	m1RingLookup[node].nodes.clear();
	m1RingLookup[node].tris.clear();
	int start = mCorners[corner].prev;
	int current = start;
	do {
		m1RingLookup[node].nodes.insert(mCorners[current].node);
		m1RingLookup[node].tris.insert(mCorners[current].tri);
		current = mCorners[mCorners[current].opposite].next;
		if (current < 0) 
			errMsg("Can't use fastNodeLookupRebuild on incomplete surfaces");
	} while (current != start);
}

void Mesh::sanityCheck(bool strict, vector<int>* deletedNodes, map<int,bool>* taintedTris) {
	const int nodes = numNodes(), tris = numTris(), corners = 3*tris;
	for(size_t i=0; i<mNodeChannels.size(); i++) {
		if (mNodeChannels[i]->size() != nodes)
			errMsg("Node channel size mismatch");
	}
	for(size_t i=0; i<mTriChannels.size(); i++) {
		if (mTriChannels[i]->size() != tris)
			errMsg("Tri channel size mismatch");
	}
	if ((int)m1RingLookup.size() != nodes)
		errMsg("1Ring size wrong");
	for(size_t t=0; t<mTris.size(); t++) { 
		if (taintedTris && taintedTris->find(t) != taintedTris->end()) continue;
		for (int c=0; c<3; c++) {
			int corner = t*3+c;
			int node = mTris[t].c[c];
			int next = mTris[t].c[(c+1)%3];
			int prev = mTris[t].c[(c+2)%3];            
			int rnext = mCorners[corner].next;
			int rprev = mCorners[corner].prev;
			int ro = mCorners[corner].opposite;
			if (node < 0 || node >= nodes || next < 0 || next >= nodes || prev < 0 || prev >= nodes)
				errMsg("invalid node entry");
			if (mCorners[corner].node != node || mCorners[corner].tri != (int)t)
				errMsg("invalid basic corner entry");
			if (rnext < 0 || rnext >= corners || rprev < 0 || rprev >= corners || ro >= corners)
				errMsg("invalid corner links");
			if (mCorners[rnext].node != next || mCorners[rprev].node != prev)
				errMsg("invalid corner next/prev");
			if (strict && ro < 0)
				errMsg("opposite missing");
			if (mCorners[ro].opposite != corner)
				errMsg("invalid opposite ref");
			set<int>& rnodes = m1RingLookup[node].nodes;
			set<int>& rtris = m1RingLookup[node].tris;
			if (rnodes.find(next) == rnodes.end() || rnodes.find(prev) == rnodes.end()) {
				debMsg("Tri "<< t << " " << node << " " << next << " " << prev , 1);
				for(set<int>::iterator it= rnodes.begin(); it != rnodes.end(); ++it)
					debMsg( *it , 1);
				errMsg("node missing in 1ring");
			}
			if (rtris.find(t) == rtris.end()) {
			   debMsg("Tri "<< t << " " << node , 1);
			   errMsg("tri missing in 1ring");
			}
		}
	}
	for (int n=0; n<nodes; n++) {
		bool docheck=true;
		if (deletedNodes)
			for (size_t e=0; e<deletedNodes->size(); e++)
				if ((*deletedNodes)[e] == n) docheck=false;;
		
		if (docheck) {
			set<int>& sn = m1RingLookup[n].nodes;
			set<int>& st = m1RingLookup[n].tris;
			set<int> sn2;
			
			for (set<int>::iterator it=st.begin(); it != st.end(); ++it) {
				bool found = false;
				for (int c=0; c<3; c++) {
					if (mTris[*it].c[c] == n)
						found = true;
					else
						sn2.insert(mTris[*it].c[c]);
				}            
				if (!found) {
					cout << *it << " " << n << endl;
					for (int c=0; c<3; c++) cout << mTris[*it].c[c] << endl;
					errMsg("invalid triangle in 1ring");
				}
				if (taintedTris && taintedTris->find(*it) != taintedTris->end()) {
					cout << *it << endl;
					errMsg("tainted tri still is use");
				}
			}
			if (sn.size() != sn2.size())
				errMsg("invalid nodes in 1ring");
			for (set<int>::iterator it=sn.begin(), it2=sn2.begin(); it != sn.end(); ++it,++it2) {
				if (*it != *it2) {
					cout << "Node " << n << ": " << *it << " vs " << *it2 << endl;
					errMsg("node ring mismatch");
				}
			}
		}        
	}
}

//*****************************************************************************
// rasterization

void meshSDF(Mesh& mesh, LevelsetGrid& levelset, Real sigma, Real cutoff = 0.);

//! helper vec3 array container
struct CVec3Ptr {
	Real *x, *y, *z; 
	inline Vec3 get(int i) const { return Vec3(x[i],y[i],z[i]); };
	inline void set(int i, const Vec3& v) { x[i]=v.x; y[i]=v.y; z[i]=v.z; };
};
//! helper vec3 array, for CUDA compatibility, remove at some point
struct CVec3Array {    
	CVec3Array(int sz) {
		x.resize(sz);
		y.resize(sz);
		z.resize(sz);        
	}    
	CVec3Array(const std::vector<Vec3>& v) {
		x.resize(v.size());
		y.resize(v.size());
		z.resize(v.size());
		for (size_t i=0; i<v.size(); i++) {
			x[i] = v[i].x;
			y[i] = v[i].y;
			z[i] = v[i].z;
		}
	}
	CVec3Ptr data() {
		CVec3Ptr a = { x.data(), y.data(), z.data()};
		return a;
	}
	inline const Vec3 operator[](int idx) const { return Vec3((Real)x[idx], (Real)y[idx], (Real)z[idx]); }
	inline void set(int idx, const Vec3& v) { x[idx] = v.x; y[idx] = v.y; z[idx] = v.z; }	
	inline int size() { return x.size(); }    
	std::vector<Real> x, y, z;
};

	
//void SDFKernel(const int* partStart, const int* partLen, CVec3Ptr pos, CVec3Ptr normal, Real* sdf, Vec3i gridRes, int intRadius, Real safeRadius2, Real cutoff2, Real isigma2);
//! helper for rasterization
static void SDFKernel(Grid<int>& partStart, Grid<int>& partLen, CVec3Ptr pos, CVec3Ptr normal, LevelsetGrid& sdf, Vec3i gridRes, int intRadius, Real safeRadius2, Real cutoff2, Real isigma2)
{
	for (int cnt_x(0); cnt_x < gridRes[0]; ++cnt_x) {
		for (int cnt_y(0); cnt_y < gridRes[1]; ++cnt_y) {
			for (int cnt_z(0); cnt_z < gridRes[2]; ++cnt_z) {
				// cell index, center
				Vec3i cell = Vec3i(cnt_x, cnt_y, cnt_z);
				if (cell.x >= gridRes.x || cell.y >= gridRes.y || cell.z >= gridRes.z) return;    
				Vec3 cpos = Vec3(cell.x + 0.5f, cell.y + 0.5f, cell.z + 0.5f);
				Real sum = 0.0f;
				Real dist = 0.0f;
				
				// query cells within block radius
				Vec3i minBlock = Vec3i(max(cell.x - intRadius,0), max(cell.y - intRadius,0), max(cell.z - intRadius,0));
				Vec3i maxBlock = Vec3i(min(cell.x + intRadius, gridRes.x - 1), min(cell.y + intRadius, gridRes.y - 1), min(cell.z + intRadius, gridRes.z - 1)); 
				for (int i=minBlock.x; i<=maxBlock.x; i++)
					for (int j=minBlock.y; j<=maxBlock.y; j++)
						for (int k=minBlock.z; k<=maxBlock.z; k++) {
							// test if block is within radius
							Vec3 d = Vec3(cell.x-i, cell.y-j, cell.z-k);                
							Real normSqr = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; 
							if (normSqr > safeRadius2) continue;
							
							// find source cell, and divide it into thread blocks
							int block = i + gridRes.x * (j + gridRes.y * k);
							int slen = partLen[block];
							if (slen == 0) continue;
							int start = partStart[block];
							
							// process sources
							for(int s=0; s<slen; s++) {                     
								
								// actual sdf kernel
								Vec3 r = cpos - pos.get(start+s);
								Real normSqr = r[0]*r[0] + r[1]*r[1] + r[2]*r[2]; 
								Real r2 = normSqr;
								if (r2 < cutoff2) {
									Real w = expf(-r2*isigma2);
									sum += w;
									dist += dot(normal.get(start+s), r) * w;                        
								}
							}
						}
				// writeback
				if (sum > 0.0f) {
					//sdf[cell.x + gridRes.x * (cell.y + gridRes.y * cell.z)] = dist / sum;    
					sdf(cell.x ,cell.y ,cell.z) = dist / sum;    
				}
			}	
		}
	}
}

static inline IndexInt _cIndex(const Vec3& pos, const Vec3i& s) {
	Vec3i p = toVec3i(pos);
	if (p.x < 0 || p.y < 0 || p.z < 0 || p.x >= s.x || p.y >= s.y || p.z >= s.z) return -1;
	return p.x + s.x * (p.y + s.y * p.z);
}

//! Kernel: Apply a shape to a grid, setting value inside
KERNEL() template<class T> 
void ApplyMeshToGrid (Grid<T>* grid, Grid<Real>& sdf, T value, FlagGrid* respectFlags) {
	if (respectFlags && respectFlags->isObstacle(i,j,k))
		return;
	if (sdf(i,j,k) < 0)
	{	
		(*grid)(i,j,k) = value;
	}
}

void Mesh::applyMeshToGrid(GridBase* grid, FlagGrid* respectFlags, Real cutoff, Real meshSigma) {
	FluidSolver dummy(grid->getSize());
	LevelsetGrid mesh_sdf(&dummy, false);
	meshSDF(*this, mesh_sdf, meshSigma, cutoff); // meshSigma=2 fixed here
	
#	if NOPYTHON!=1
	if (grid->getType() & GridBase::TypeInt)
		ApplyMeshToGrid<int> ((Grid<int>*)grid, mesh_sdf, _args.get<int>("value"), respectFlags);
	else if (grid->getType() & GridBase::TypeReal)
		ApplyMeshToGrid<Real> ((Grid<Real>*)grid, mesh_sdf, _args.get<Real>("value"), respectFlags);
	else if (grid->getType() & GridBase::TypeVec3)
		ApplyMeshToGrid<Vec3> ((Grid<Vec3>*)grid, mesh_sdf, _args.get<Vec3>("value"), respectFlags);
	else
		errMsg("Shape::applyToGrid(): unknown grid type");
#	else
	errMsg("Not yet supported...");
#	endif
}

void Mesh::computeLevelset(LevelsetGrid& levelset, Real sigma, Real cutoff) {
	meshSDF( *this, levelset, sigma, cutoff); 
}

LevelsetGrid Mesh::getLevelset(Real sigma, Real cutoff) {
	LevelsetGrid phi(getParent());
	meshSDF(*this, phi, sigma, cutoff);
	return phi;
}

void meshSDF(Mesh& mesh, LevelsetGrid& levelset, Real sigma, Real cutoff)
{  
	if (cutoff<0) cutoff = 2*sigma;
	Real maxEdgeLength = 0.75;
	Real numSamplesPerCell = 0.75;
	
	Vec3i gridRes = levelset.getSize();
	Vec3 mult = toVec3(gridRes) / toVec3(mesh.getParent()->getGridSize());
	
	// prepare center values
	std::vector<Vec3> center;
	std::vector<Vec3> normals;
	short bigEdges(0);
	std::vector<Vec3> samplePoints;
	for(int i=0; i<mesh.numTris(); i++){	
		center.push_back(Vec3(mesh.getFaceCenter(i) * mult));
		normals.push_back(mesh.getFaceNormal(i));
		//count big, stretched edges
		bigEdges = 0;
		for (short edge(0); edge <3; ++edge){
			if(norm(mesh.getEdge(i,edge)) > maxEdgeLength){
				bigEdges += 1 << edge;
			}
		}
		if(bigEdges > 0){
			samplePoints.clear();
			short iterA, pointA, iterB, pointB;
			int numSamples0 = norm(mesh.getEdge(i,1)) * numSamplesPerCell;
			int numSamples1 = norm(mesh.getEdge(i,2)) * numSamplesPerCell;
			int numSamples2 = norm(mesh.getEdge(i,0)) * numSamplesPerCell;
			if(! (bigEdges & (1 << 0))){
				//loop through 0,1
				iterA = numSamples1;
				pointA = 0;
				iterB = numSamples2;
				pointB = 1;
			}
			else if(! (bigEdges & (1 << 1))){
				//loop through 1,2
				iterA = numSamples2;
				pointA = 1;
				iterB = numSamples0;
				pointB = 2;
				
			}
			else{
				//loop through 2,0
				iterA = numSamples0;
				pointA = 2;
				iterB = numSamples1;
				pointB = 0;
			}
			
			Real u(0.),v(0.),w(0.); // barycentric uvw coords
			Vec3 samplePoint,normal;
			for (int sample0(0); sample0 < iterA; ++sample0){
				u = Real(1. * sample0 / iterA);
				for (int sample1(0); sample1 < iterB; ++sample1){
					v = Real(1. * sample1 / iterB);
					w = 1 - u - v;
					if (w < 0.)
						continue;
					samplePoint = 	mesh.getNode(i,pointA) * mult * u + 
					mesh.getNode(i,pointB) * mult * v + 
					mesh.getNode(i,(3 - pointA - pointB)) * mult * w;
					samplePoints.push_back(samplePoint);
					normal = mesh.getFaceNormal(i);
					normals.push_back(normal);							
				}
			}
			center.insert(center.end(), samplePoints.begin(), samplePoints.end());
		}
	}	

	// prepare grid    
	levelset.setConst( -cutoff );
	
	// 1. count sources per cell
	Grid<int> srcPerCell(levelset.getParent());
	for (size_t i=0; i<center.size(); i++) {
		IndexInt idx = _cIndex(center[i], gridRes);
		if (idx >= 0)
			srcPerCell[idx]++;
	}
	
	// 2. create start index lookup
	Grid<int> srcCellStart(levelset.getParent());
	int cnt=0;
	FOR_IJK(srcCellStart) {
		IndexInt idx = srcCellStart.index(i,j,k);
		srcCellStart[idx] = cnt;
		cnt += srcPerCell[idx];
	}
	
	// 3. reorder nodes
	CVec3Array reorderPos(center.size());
	CVec3Array reorderNormal(center.size());
	{
		Grid<int> curSrcCell(levelset.getParent());
		for (int i=0; i<(int)center.size(); i++) {
			IndexInt idx = _cIndex(center[i], gridRes);
			if (idx < 0) continue;
			IndexInt idx2 = srcCellStart[idx] + curSrcCell[idx];
			reorderPos.set(idx2, center[i]);
			reorderNormal.set(idx2, normals[i]);
			curSrcCell[idx]++;
		}
	}
	
	// construct parameters
	Real safeRadius = cutoff + sqrt(3.0)*0.5;
	Real safeRadius2 = safeRadius*safeRadius;
	Real cutoff2 = cutoff*cutoff;
	Real isigma2 = 1.0/(sigma*sigma);
	int intRadius = (int)(cutoff+0.5);
	
	SDFKernel( srcCellStart, srcPerCell,
			  reorderPos.data(), reorderNormal.data(), 
			  levelset, gridRes, intRadius, safeRadius2, cutoff2, isigma2);
	
	// floodfill outside
	std::stack<Vec3i> outside;
	FOR_IJK(levelset) {
		if (levelset(i,j,k) >= cutoff-1.0f) 
			outside.push(Vec3i(i,j,k));
	}
	while(!outside.empty()) {
		Vec3i c = outside.top();
		outside.pop();
		levelset(c) = cutoff;
		if (c.x > 0 && levelset(c.x-1, c.y, c.z) < 0) outside.push(Vec3i(c.x-1,c.y,c.z));
		if (c.y > 0 && levelset(c.x, c.y-1, c.z) < 0) outside.push(Vec3i(c.x,c.y-1,c.z));
		if (c.z > 0 && levelset(c.x, c.y, c.z-1) < 0) outside.push(Vec3i(c.x,c.y,c.z-1));
		if (c.x < levelset.getSizeX()-1 && levelset(c.x+1, c.y, c.z) < 0) outside.push(Vec3i(c.x+1,c.y,c.z));
		if (c.y < levelset.getSizeY()-1 && levelset(c.x, c.y+1, c.z) < 0) outside.push(Vec3i(c.x,c.y+1,c.z));
		if (c.z < levelset.getSizeZ()-1 && levelset(c.x, c.y, c.z+1) < 0) outside.push(Vec3i(c.x,c.y,c.z+1));
	};
}
	
// Blender data pointer accessors
std::string Mesh::getNodesDataPointer() {
	std::ostringstream out;
	out << &mNodes;
	return out.str();
}
std::string Mesh::getTrisDataPointer() {
	std::ostringstream out;
	out << &mTris;
	return out.str();
}

// mesh data

MeshDataBase::MeshDataBase(FluidSolver* parent) :
		PbClass(parent) , mMesh(NULL) {
}

MeshDataBase::~MeshDataBase()
{
	// notify parent of deletion
	if(mMesh)
		mMesh->deregister(this);
}

// actual data implementation

template<class T>
MeshDataImpl<T>::MeshDataImpl(FluidSolver* parent) :
	MeshDataBase(parent) , mpGridSource(NULL), mGridSourceMAC(false) {
}

template<class T>
MeshDataImpl<T>::MeshDataImpl(FluidSolver* parent, MeshDataImpl<T>* other) :
	MeshDataBase(parent) , mpGridSource(NULL), mGridSourceMAC(false) {
	this->mData = other->mData;
	setName(other->getName());
}

template<class T>
MeshDataImpl<T>::~MeshDataImpl() {
}

template<class T>
IndexInt MeshDataImpl<T>::getSizeSlow() const {
	return mData.size();
}
template<class T>
void MeshDataImpl<T>::addEntry() {
	// add zero'ed entry
	T tmp = T(0.);
	// for debugging, force init:
	//tmp = T(0.02 * mData.size()); // increasing
	//tmp = T(1.); // constant 1
	return mData.push_back(tmp);
}
template<class T>
void MeshDataImpl<T>::resize(IndexInt s) {
	mData.resize(s);
}
template<class T>
void MeshDataImpl<T>::copyValueSlow(IndexInt from, IndexInt to) {
	this->copyValue(from,to);
}
template<class T>
MeshDataBase* MeshDataImpl<T>::clone() {
	MeshDataImpl<T>* npd = new MeshDataImpl<T>( getParent(), this );
	return npd;
}

template<class T>
void MeshDataImpl<T>::setSource(Grid<T>* grid, bool isMAC ) {
	mpGridSource = grid;
	mGridSourceMAC = isMAC;
	if(isMAC) assertMsg( dynamic_cast<MACGrid*>(grid) != NULL , "Given grid is not a valid MAC grid");
}

template<class T>
void MeshDataImpl<T>::initNewValue(IndexInt idx, Vec3 pos) {
	if(!mpGridSource)
		mData[idx] = 0;
	else {
		mData[idx] = mpGridSource->getInterpolated(pos);
	}
}

// special handling needed for velocities
template<>
void MeshDataImpl<Vec3>::initNewValue(IndexInt idx, Vec3 pos) {
	if(!mpGridSource)
		mData[idx] = 0;
	else {
		if(!mGridSourceMAC)
			mData[idx] = mpGridSource->getInterpolated(pos);
		else
			mData[idx] = ((MACGrid*)mpGridSource)->getInterpolated(pos);
	}
}

//! update additional mesh data
void Mesh::updateDataFields() {
	for(size_t i=0; i<mNodes.size(); ++i) {
		Vec3 pos = mNodes[i].pos;
		for(IndexInt md=0; md<(IndexInt)mMdataReal.size(); ++md)
			mMdataReal[md]->initNewValue(i, mNodes[i].pos);
		for(IndexInt md=0; md<(IndexInt)mMdataVec3.size(); ++md)
			mMdataVec3[md]->initNewValue(i, mNodes[i].pos);
		for(IndexInt md=0; md<(IndexInt)mMdataInt.size(); ++md)
			mMdataInt[md]->initNewValue(i, mNodes[i].pos);
	}
}

template<typename T>
void MeshDataImpl<T>::load(string name) {
	if (name.find_last_of('.') == string::npos)
		errMsg("file '" + name + "' does not have an extension");
	string ext = name.substr(name.find_last_of('.'));
	if ( ext == ".uni")
		readMdataUni<T>(name, this);
	else if ( ext == ".raw") // raw = uni for now
		readMdataUni<T>(name, this);
	else
		errMsg("mesh data '" + name +"' filetype not supported for loading");
}

template<typename T>
void MeshDataImpl<T>::save(string name) {
	if (name.find_last_of('.') == string::npos)
		errMsg("file '" + name + "' does not have an extension");
	string ext = name.substr(name.find_last_of('.'));
	if (ext == ".uni")
		writeMdataUni<T>(name, this);
	else if (ext == ".raw") // raw = uni for now
		writeMdataUni<T>(name, this);
	else
		errMsg("mesh data '" + name +"' filetype not supported for saving");
}

// specializations

template<>
MeshDataBase::MdataType MeshDataImpl<Real>::getType() const {
	return MeshDataBase::TypeReal;
}
template<>
MeshDataBase::MdataType MeshDataImpl<int>::getType() const {
	return MeshDataBase::TypeInt;
}
template<>
MeshDataBase::MdataType MeshDataImpl<Vec3>::getType() const {
	return MeshDataBase::TypeVec3;
}

KERNEL(pts) template<class T> void knSetMdataConst(MeshDataImpl<T>& mdata, T value) { mdata[idx] = value; }

KERNEL(pts) template<class T, class S> void knMdataSet  (MeshDataImpl<T>& me, const MeshDataImpl<S>& other) { me[idx] += other[idx]; }
KERNEL(pts) template<class T, class S> void knMdataAdd  (MeshDataImpl<T>& me, const MeshDataImpl<S>& other) { me[idx] += other[idx]; }
KERNEL(pts) template<class T, class S> void knMdataSub  (MeshDataImpl<T>& me, const MeshDataImpl<S>& other) { me[idx] -= other[idx]; }
KERNEL(pts) template<class T, class S> void knMdataMult (MeshDataImpl<T>& me, const MeshDataImpl<S>& other) { me[idx] *= other[idx]; }
KERNEL(pts) template<class T, class S> void knMdataDiv  (MeshDataImpl<T>& me, const MeshDataImpl<S>& other) { me[idx] /= other[idx]; }

KERNEL(pts) template<class T, class S> void knMdataSetScalar (MeshDataImpl<T>& me, const S& other)  { me[idx]  = other; }
KERNEL(pts) template<class T, class S> void knMdataAddScalar (MeshDataImpl<T>& me, const S& other)  { me[idx] += other; }
KERNEL(pts) template<class T, class S> void knMdataMultScalar(MeshDataImpl<T>& me, const S& other)  { me[idx] *= other; }
KERNEL(pts) template<class T, class S> void knMdataScaledAdd (MeshDataImpl<T>& me, const MeshDataImpl<T>& other, const S& factor) { me[idx] += factor * other[idx]; }

KERNEL(pts) template<class T> void knMdataSafeDiv (MeshDataImpl<T>& me, const MeshDataImpl<T>& other) { me[idx] = safeDivide(me[idx], other[idx]); }
KERNEL(pts) template<class T> void knMdataSetConst(MeshDataImpl<T>& mdata, T value) { mdata[idx] = value; }

KERNEL(pts) template<class T> void knMdataClamp (MeshDataImpl<T>& me, T min, T max) { me[idx] = clamp( me[idx], min, max); }
KERNEL(pts) template<class T> void knMdataClampMin(MeshDataImpl<T>& me, const T vmin)               { me[idx] = std::max(vmin, me[idx]); }
KERNEL(pts) template<class T> void knMdataClampMax(MeshDataImpl<T>& me, const T vmax)               { me[idx] = std::min(vmax, me[idx]); }
KERNEL(pts)                   void knMdataClampMinVec3(MeshDataImpl<Vec3>& me, const Real vmin)     {
	me[idx].x = std::max(vmin, me[idx].x);
	me[idx].y = std::max(vmin, me[idx].y);
	me[idx].z = std::max(vmin, me[idx].z);
}
KERNEL(pts)                   void knMdataClampMaxVec3(MeshDataImpl<Vec3>& me, const Real vmax)     {
	me[idx].x = std::min(vmax, me[idx].x);
	me[idx].y = std::min(vmax, me[idx].y);
	me[idx].z = std::min(vmax, me[idx].z);
}

// python operators


template<typename T>
MeshDataImpl<T>& MeshDataImpl<T>::copyFrom(const MeshDataImpl<T>& a) {
	assertMsg (a.mData.size() == mData.size() , "different mdata size "<<a.mData.size()<<" vs "<<this->mData.size() );
	memcpy( &mData[0], &a.mData[0], sizeof(T) * mData.size() );
	return *this;
}

template<typename T>
void MeshDataImpl<T>::setConst(T s) {
	knMdataSetScalar<T,T> op( *this, s );
}

template<typename T>
void MeshDataImpl<T>::setConstRange(T s, const int begin, const int end) {
	for(int i=begin; i<end; ++i) (*this)[i] = s;
}

// special set by flag
KERNEL(pts) template<class T, class S> void knMdataSetScalarIntFlag(MeshDataImpl<T>& me, const S& other, const MeshDataImpl<int>& t, const int itype) {
	if(t[idx]&itype) me[idx] = other;
}
template<typename T>
void MeshDataImpl<T>::setConstIntFlag(T s, const MeshDataImpl<int>& t, const int itype) {
	knMdataSetScalarIntFlag<T,T> op(*this, s, t, itype);
}

template<typename T>
void MeshDataImpl<T>::add(const MeshDataImpl<T>& a) {
	knMdataAdd<T,T> op( *this, a );
}
template<typename T>
void MeshDataImpl<T>::sub(const MeshDataImpl<T>& a) {
	knMdataSub<T,T> op( *this, a );
}

template<typename T>
void MeshDataImpl<T>::addConst(T s) {
	knMdataAddScalar<T,T> op( *this, s );
}

template<typename T>
void MeshDataImpl<T>::addScaled(const MeshDataImpl<T>& a, const T& factor) {
	knMdataScaledAdd<T,T> op( *this, a, factor );
}

template<typename T>
void MeshDataImpl<T>::mult( const MeshDataImpl<T>& a) {
	knMdataMult<T,T> op( *this, a );
}

template<typename T>
void MeshDataImpl<T>::safeDiv(const MeshDataImpl<T>& a) {
	knMdataSafeDiv<T> op( *this, a );
}

template<typename T>
void MeshDataImpl<T>::multConst(T s) {
	knMdataMultScalar<T,T> op( *this, s );
}


template<typename T>
void MeshDataImpl<T>::clamp(Real vmin, Real vmax) {
	knMdataClamp<T> op( *this, vmin, vmax );
}

template<typename T>
void MeshDataImpl<T>::clampMin(Real vmin) {
	knMdataClampMin<T> op( *this, vmin );
}
template<typename T>
void MeshDataImpl<T>::clampMax(Real vmax) {
	knMdataClampMax<T> op( *this, vmax );
}

template<>
void MeshDataImpl<Vec3>::clampMin(Real vmin) {
	knMdataClampMinVec3 op( *this, vmin );
}
template<>
void MeshDataImpl<Vec3>::clampMax(Real vmax) {
	knMdataClampMaxVec3 op( *this, vmax );
}

template<typename T> KERNEL(pts, reduce=+) returns(T result=T(0.)) T    KnPtsSum(const MeshDataImpl<T>& val, const MeshDataImpl<int> *t, const int itype) { if(t && !((*t)[idx]&itype)) return; result += val[idx]; }
template<typename T> KERNEL(pts, reduce=+) returns(Real result=0.) Real KnPtsSumSquare(const MeshDataImpl<T>& val)    { result += normSquare(val[idx]); }
template<typename T> KERNEL(pts, reduce=+) returns(Real result=0.) Real KnPtsSumMagnitude(const MeshDataImpl<T>& val) { result += norm(val[idx]); }

template<typename T>
T MeshDataImpl<T>::sum(const MeshDataImpl<int> *t, const int itype) const {
	return KnPtsSum<T>(*this, t, itype);
}
template<typename T>
Real MeshDataImpl<T>::sumSquare() const {
	return KnPtsSumSquare<T>(*this);
}
template<typename T>
Real MeshDataImpl<T>::sumMagnitude() const {
	return KnPtsSumMagnitude<T>(*this);
}

template<typename T>
KERNEL(pts, reduce=min) returns(Real minVal=std::numeric_limits<Real>::max())
Real CompMdata_Min(const MeshDataImpl<T>& val) {
	if (val[idx] < minVal)
		minVal = val[idx];
}

template<typename T>
KERNEL(pts, reduce=max) returns(Real maxVal=-std::numeric_limits<Real>::max())
Real CompMdata_Max(const MeshDataImpl<T>& val) {
	if (val[idx] > maxVal)
		maxVal = val[idx];
}

template<typename T>
Real MeshDataImpl<T>::getMin() {
	return CompMdata_Min<T> (*this);
}

template<typename T>
Real MeshDataImpl<T>::getMaxAbs() {
	Real amin = CompMdata_Min<T> (*this);
	Real amax = CompMdata_Max<T> (*this);
	return max( fabs(amin), fabs(amax));
}

template<typename T>
Real MeshDataImpl<T>::getMax() {
	return CompMdata_Max<T> (*this);
}

template<typename T>
void MeshDataImpl<T>::printMdata(IndexInt start, IndexInt stop, bool printIndex)
{
	std::ostringstream sstr;
	IndexInt s = (start>0 ? start : 0                      );
	IndexInt e = (stop>0  ? stop  : (IndexInt)mData.size() );
	s = Manta::clamp(s, (IndexInt)0, (IndexInt)mData.size());
	e = Manta::clamp(e, (IndexInt)0, (IndexInt)mData.size());

	for(IndexInt i=s; i<e; ++i) {
		if(printIndex) sstr << i<<": ";
		sstr<<mData[i]<<" "<<"\n";
	}
	debMsg( sstr.str() , 1 );
}
template<class T> std::string MeshDataImpl<T>::getDataPointer() {
	std::ostringstream out;
	out << &mData;
	return out.str();
}

// specials for vec3
// work on length values, ie, always positive (in contrast to scalar versions above)

KERNEL(pts, reduce=min) returns(Real minVal=-std::numeric_limits<Real>::max())
Real CompMdata_MinVec3(const MeshDataImpl<Vec3>& val) {
	const Real s = normSquare(val[idx]);
	if (s < minVal)
		minVal = s;
}

KERNEL(pts, reduce=max) returns(Real maxVal=-std::numeric_limits<Real>::min())
Real CompMdata_MaxVec3(const MeshDataImpl<Vec3>& val) {
	const Real s = normSquare(val[idx]);
	if (s > maxVal)
		maxVal = s;
}

template<>
Real MeshDataImpl<Vec3>::getMin() {
	return sqrt(CompMdata_MinVec3 (*this));
}

template<>
Real MeshDataImpl<Vec3>::getMaxAbs() {
	return sqrt(CompMdata_MaxVec3 (*this));  // no minimum necessary here
}

template<>
Real MeshDataImpl<Vec3>::getMax() {
	return sqrt(CompMdata_MaxVec3 (*this));
}


// explicit instantiation
template class MeshDataImpl<int>;
template class MeshDataImpl<Real>;
template class MeshDataImpl<Vec3>;


} //namespace
